Fictitious time wave packet dynamics: II. Hydrogen atom in external fields 
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In the preceding paper [T. Fabcic et al., preprint] "restricted Gaussian wave packets" were in- 
troduced for the regularized Coulomb problem in the four-dimensional Kustaanheimo-Stiefel coor- 
dinates, and their exact time propagation was derived analytically in a fictitious time variable. We 
now establish the Gaussian wave packet method for the hydrogen atom in static external fields. 
A superposition of restricted Gaussian wave packets is used as a trial function in the application 
of the time-dependent variational principle. The external fields introduce couplings between the 
basis states. The set of coupled wave packets is propagated numerically, and eigenvalues of the 
Schrodinger equation are obtained by the frequency analysis of the time autocorrelation function. 
The advantage of the wave packet propagation in the fictitious time variable is that the computations 
are exact for the field-free hydrogen atom and approximations from the time-dependent variational 
principle only stem from the external fields. Examples are presented for the hydrogen atom in a 
magnetic field and in crossed electric and magnetic fields. 

PACS numbers: 32.80.Ee, 31.15.xt, 32.60.+i, 05.45.-a 
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I. INTRODUCTION 

The hydrogen atom in a static magnetic field [J, [H, 
and in crossed electric and magnetic fields 
I M, El E H EI El 13 is a non-integrable 
system which can be accessed both experimentally and 
theoretically and has attracted much attention during re- 
cent decades. Exact quantum spectra of the system can 
be obtained by numerical diagonalization of the Hamil- 
tonian in a large Sturmian type basis set. Neverthe- 
less, the atom has served as an example system for the 
development and verification of alternative quantization 
methods, e.g., semiclassical closed-orbit theory [l7l fl8j. 
periodic-orbit theory [H [2(| , and cycle-expansion tech- 
niques [2l| . 

Another alternative to large quantum computations is 
the application of the time-dependent variational princi- 
ple (TDVP) [22|]. For a wave packet depending on a set of 
variational parameters the timc-dcpcndcnt Schrodinger 
equation is transformed to a system of ordinary differen- 
tial equations for the variational parameters. Quantum 
spectra can be obtained by a frequency analysis of the 
time autocorrelation function of the wave packet. The 
method has been established by Heller [H [24| for single 
or coupled Gaussian wave packets (GWPs). It is well 
suited for nonsingular smooth potentials but certainly 
far from ideal for atomic systems with singular Coulomb 
potentials. 

The wave packet dynamics in atomic systems has been 
studied for the field-free hydrogen atom [H [H \2l\ . 
and in particular for the atom in time-dependent exter- 
nal fields, e.g., microwaves or short laser pulses. While 
Rydberg wave packets are usually dispersive, the possi- 
ble existence of nondispersive coherent states has been 
demonstrated for the hydrogen atom in microwave fields 

MM- 

In the preceding paper [3(| we have established the 
Gaussian wave packet method for the Coulomb problem. 



Using the Kustaanheimo-Stiefel (KS) regularization the 
singular Coulomb problem was transformed to the four- 
dimensional (4D) harmonic oscillator with a constraint. 
We introduced the set of "restricted Gaussian wave pack- 
ets" obeying that constraint by confining the space of 
the Gaussian parameters. The exact propagation of the 
restricted GWPs in a fictitious time variable could be 
derived analytically. 

In this paper we extend the fictitious time wave packet 
propagation to the hydrogen atom in static external elec- 
tric and magnetic fields. A superposition of restricted 
GWPs is used as the variational trial function. The 
time-dependent variational principle is applied in such a 
way that the wave packet dynamics is exact for the field- 
free hydrogen atom and couplings between the GWPs are 
only induced by the external fields. In the presence of a 
single external homogeneous field the rotational symme- 
try of the hydrogen atom is preserved and one component 
of the angular momentum, say l z , is conserved. In that 
case we employ the modified 2D Gaussian wave pack- 
ets with well-defined magnetic quantum number m intro- 
duced and discussed in Ref. [30( and perform computa- 
tions in the subspaces of the different magnetic quantum 
numbers m separately. In crossed fields the cylindrical 
symmetry is broken and computations are performed in 
the basis of the restricted 3D GWPs without well-defined 
angular momentum quantum numbers. 

The fact that the wave packet propagation is exact for 
the pure Coulomb problem might imply that the external 
fields are treated as a perturbation and the method does 
not work well beyond the perturbative regime. However, 
this is not the case. The dynamics of wave packets is 
exact to all orders in the field strengths within the al- 
lowed set of trial wave functions, i.e., the variational ap- 
proximation only concerns the restriction of the Hilbcrt 
space. The power of the method will be demonstrated 
by application to the diamagnetic hydrogen atom in the 
strong non-pertubative regime at the field-free ionization 
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threshold. 

The paper is organized as follows. In Sec. [IT] we intro- 
duce the regularization and scaling of the Hamiltonian 
with external fields and discuss the general idea of how 
to obtain quantum spectra by frequency analysis of the 
fictitious time autocorrelation function of the propagated 
wave packets. In Sec. IHII the time-dependent variational 
principle is explained. The equations of motion for the 
variational parameters are derived for the superposition 
of restricted 3D and modified 2D GWPs, and the numeri- 
cal time propagation of coupled wave packets is discussed. 
Results for the diamagnetic hydrogen atom and the atom 
in crossed electric and magnetic fields are presented in 
Sec. |IV] Concluding remarks are given in Sec. |V] 



II. REGULARIZED HYDROGEN ATOM IN 
EXTERNAL FIELDS 

In the preceding paper [30( the fictitious time wave 
packet dynamics has been discussed for the field-free hy- 
drogen atom. We now consider the atom in external elec- 
tric and magnetic fields. For perpendicular fields with the 
electric and magnetic field along the x and z axis, respec- 
tively, the Hamiltonian in the three-dimensional coordi- 
nates reads (in atomic units with F$ = 5.14 x 10 9 V/cm, 
B Q = 2.35 x 10 5 T) 



ft - \v> 



1 1 



1 



Bl z + -B 2 {x 2 + y 2 ) + Fx . (1) 



The starting point for our investigations is the 
Schrodinger equation in the 4D Kustaanheimo-Stiefel co- 
ordinates u with x = U1U3 — 1*21*4, y = 1*11*4 + 1*21*3, and 

„,2 „.2\ T„. 

2V^1 1 ^2 

and momenta u 
the procedure of Sec. II in |3C 



2 = 5(1*1 +u\ — iff — 1*4). Introducing scaled coordinates 



1/2 



-1/2 



"cff 



p u and following 



we obtain 



1 2 
2 P " 



1 



n 2 cS E + -{n 2 cS B) 2 (u\ + u 2 2 )(u 2 3 + u\) 

l s F(mu 3 - U2U4)] u 2 



-7* 2 ff B [(Uxp 2 - U 2 pi)(ul + u\) 

+ (u 3 p4 - u 4 p 3 )(uj + I* 2 )] \ip = 2n eS ip . 



(2) 



In KS coordinates physical wave functions must fulfill the 
constraint 

(i* 2 pi - uip 2 - u 4 p 3 + U3P4) ip = . (3) 

By choosing constant parameters 

a = -n 2 cB E, P = n 2 s B, ( = n 3 cS F , (4) 

Eq. ([2]) becomes an eigenvalue problem for the effective 
quantum number n c fj . For a set of parameters (a, (3, C) 
and a given eigenvalue n e ff the energy and field strengths 
of the physical state are obtained from Eq. Q. The 



quantized energies and field strengths are located on lines 
with constant E/B and E/F 2 / 3 . 

In analogy with the field-free hydrogen atom in [3(| we 
can extend Eq. @ to the time-dependent Schrodinger 
equation in the dimensionlcss fictitious time t by the 
replacement 2n e ff — ► ijp, viz. 



i> = ( ^vl + v ) v = ( -Ja 4 + v ) i> = (t + v)i> , 



8 
'd-r' 

where V is defined via Eq. ^ as the sum of a harmonic 
potential and the contributions of the external fields. For 
the field- free hydrogen atom, i.e., a harmonic potential 
V, wave packets can be propagated analytically in the 
fictitious time [30j . 

Our goal for the hydrogen atom in external fields is 
to compute the propagation of an initial wave packet 
V'(O) by applying the time-dependent variational princi- 
ple. To this end the wave function is assumed to depend 
on a set of appropriately chosen parameters whose time- 
dependence is obtained by solving ordinary differential 
equations. The ansatz for the wave function depends on 
the symmetry of the problem. For the hydrogen atom in 
crossed fields we choose a superposition of iV restricted 
Gaussian wave packets [30[ 



V(r) 



N 

E< 

k=l 



( i[uA fc (T)u+7 fe (r)] 



(6) 



with the symmetric width matrices 



A = 









ft,, 







(7) 



y 



-a, 



depending on the four parameters (a^^a^^ax^y), and 
with 7 determining the normalization and phase of the re- 
stricted GWPs. The special form of the ansatz ©, which 
depends on, in total, 5iV time-dependent variational pa- 
rameters (instead of 15 AT complex parameters for the 
most general superposition of Gaussian wave packets in a 
4D coordinate space) guarantees that the wave function 
obeys the constraint ©. 

The hydrogen atom in a pure magnetic field, i.e., £ = 0, 
is cylindrically symmetric around the z axis, and the an- 
gular momentum component l z — m is an exact quantum 
number. For wave packets with given m quantum num- 
ber we use the ansatz 



iA k (T)M+f k (T)} e imtp 



^ m (r) = (/u/)M£y^ 
fc=i 

N 

= (pv)\ m \ e i [ a f''(' r )^ 2+a ^ r ) I/2 + 7fe ( r )]e imv ,(8) 



k=l 



with the diagonal form of the matrix A obtained by set- 
ting a x = a v = in |(7J), and semiparabolic coordinates 
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/i = \/u\ + u\ = \Jr + z and v = \Ju\ + u 
are introduced. The wave function © thus depends on 
a set of SN time-dependent variational parameters. As 
the paramagnetic term is constant this term can be ab- 
sorbed by an energy shift E — > E' = E — mB/2. In 
scmiparabolic coordinates the kinetic and potential term 
in ([5]) for the diamagnetic hydrogen atom then take the 
form 



T 
V 



ld_ 
1 



dv 2 



ld_ 

v dv 



(9) 



Once the time-dependent wave packets © or © 

are determined the eigenvalues n^2 of the stationary 
Schrodinger equation (|2|) and thus quantum spectra of 
the hydrogen atom in external fields are obtained by a 
frequency analysis of the time signal 



C(r) = W(OMt)) =5>je 



-V2n~ 



(10) 



with the amplitudes Cj depending on the choice of the 
initial wave packet. The advantage of using the ficti- 
tious time r is that the computations are exact for the 
field-free hydrogen atom and approximations from the 
time-dependent variational principle only stem from the 
external fields. By contrast, wave packet propagation in 
the physical time t is a very nontrivial task even for the 
pure Coulomb potential. 



III. TIME-DEPENDENT VARIATIONAL 
PRINCIPLE 

The propagation of the wave packets investigated in 
this paper is based on the application of the time- 
dependent variational principle. For the convenience of 
the reader we first give a brief general introduction to the 
TDVP which is then applied to the special form of the 
trial functions ([6]) and ([8]) . The formulation of McLachlan 
[131, or equivalently the minimum error method [3lj |. re- 
quires the norm of the deviation between the right-hand 
and the left-hand side of the time-dependent Schrodinger 
equation to be minimized with respect to the trial func- 
tion. The quantity 



I = \\i<p{t) - Hi){t)\\ 2 = min 



(11) 



is to be varied with respect to <p only, and then ip = <p 
is chosen, i.e., for any time t the fixed wave function 
ip(t) is supposed to be given and its time derivative ip{t) 
is determined by the requirement to minimize I. The 
equality I = is provided by the exact solution of the 
Schrodinger equation, while / in general takes positive 
values if ip is constrained by the functional form of ip. 
The wave function ip(t) is assumed to be parametrized 
by a set of complex parameters z(t) = (zi(t), . . . , z np (t)), 




FIG. 1: Sketch of the manifold M of approximation of the 
trial wave function tp(z). The variational evolution of the trial 
function, denoted by the arrow with the white arrowhead, is 
obtained as the projection of the exact time evolution —iHip, 
denoted by the arrow with the black arrowhead, onto the 
tangent space T^M of the manifold M in the point ip. 



ip(t) = ip(z(t)). For brevity the arguments of the wave 
function are dropped in the following. For parametrized 
trial functions the variations 8<p carry over to variations 
Si and the variation leads to the equations of motion 



iip - Hip ) = , 



dip 

which can be written in matrix form 



(12) 



Kz 



-ih with K = ( —— 



dtp 

l)z~ 



Hip 



dip\ /dip 
dz I \dz 

(13) 

An illustration of Eq. (| 1 2|) is presented in Fig. [TJ Here 
the manifold of approximation M , consisting of all pos- 
sible configurations ip(z), is plotted schematically as a 
2D-surface in the Hilbert space. The tangent space of 
the manifold in the point ip is a linear vector space and 
is spanned by the derivatives J-^-, k = l,...,n p . The 
tangent space is denoted by T^M in Fig. [TJ According 
to the Schrodinger equation the exact time derivative ip is 
given by —iHtp, denoted by the arrow with the black ar- 
rowhead. In general the exact time derivative does not lie 
in the tangent space, otherwise the trial function would 
be an exact solution of the Schrodinger equation. The 
variational approximation to the exact time derivative is 
given by that vector of the tangent space which has mini- 
mal deviation from the exact one. This is the orthogonal 
projection of the exact time derivative onto the tangent 
space, denoted by the arrow with the white arrowhead in 

Fig.m 

For parametrized wave functions the variational prin- 
ciple (jll[) simply reduces to a quadratic minimization 
problem where the gradient of I with respect to the time 
derivatives of the parameters must be zero 



dl 

dz k 







1. 



. , Up 



(14) 



and the TDVP leads to a reduction of the Schrodinger 
equation to a system of ordinary first-order differential 
equations of motion for the parameters z{t). The matrix 
equation (|13|) must be solved numerically after each time 
step of integration for the time derivatives z if a numerical 
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algorithm for ordinary differential equations, e.g. Runge- 
Kutta or Adams, is used. 

We now apply the time-dependent variational princi- 
ple, first in Sec. If 11 Al to the trial function ||5J) of the dia- 
magnetic hydrogen atom, and then in Sec. IIII Bl to the 
trial function (J6j> of the hydrogen atom in crossed electric 
and magnetic fields. For Gaussian type trial functions it 
is convenient to split the Hamiltonian into the kinetic 
and potential part, i.e., H = T + V , and to apply Eq. 
dHJ) in the form 



dip 



iip — Tip 



dtp 
dz 



Vip 



(15) 



Note that the variational approach substantially differs 
from a perturbative treatment of the hydrogen atom 
in external fields, and is valid even in the strong non- 
perturbative regime. 



A. Diamagnetic hydrogen atom 

For the time-dependent wave packets of the hydro- 
gen atom in a homogeneous external magnetic field with 
given m quantum number we use the ansatz (jHJ) which 
can be written in the form 



N 



with the basis states 

9m (y) = 0H |r 



fe=i 



J-ia^n +a,vv +7) 



(16) 



(17) 



As already mentioned, the cylindrical symmetry of the 
system is accounted for by setting a x — a y = and only 



the time-dependent parameters z = (y , 



with 



y = (7, a M , a„) remain. The evolution of the basis states 
is obtained by the TDVP. The variational equations of 
motion are set up by evaluating Eq. (fl~5)) . First we let the 
time derivative and the Laplacian act on the basis states 
(fTf) to obtain 

[-j k + 2i (o*+o*)(l + H)- 
(d* + 2(a*) V - (d k u + 2(at) 2 y] 5m (y fc ,x) 



5m(y fe ,x) 



(18) 



for k = 1, . . . , N. Eq. JTS]) defines the coefficients v k , V k , 
V k as functions of the parameters a k , a k and the time 
derivatives 7 , a k , a k . The equations of motion can be 
written as 



-2{atf 



-V k 

2 " ' 



1 . 



-V* 



(19a) 
(19b) 



j k = 2i(a k +at)(l + \m\)-v k , (19c) 
1, . . . , TV, and the yet unknown coefficients V k , 
v k . Note that the equations of motion (flU|) 



with k = 
V k , and 

are in general coupled through the coefficients vfi, ... 
V k which become time-dependent in the presence of an- 
harmonic potentials. They must be determined from a 
system of linear equations, which follows from Eq. (|15|) 
when inserting the trial function (|16[) . Using the deriva- 
tives of the basis states (|17l) with respect to the varia- 
tional parameters, viz. -^g k n = ig k ml -£rrgt = *M 2 ffm» 



and ^g k , 



mi E q- 

the matrix equation 



U) of the TDVP finally yields 



N 



E WJgtH + \{g l m \Ag m )vH + \{g l m \Ag m )v k 

k=l ^ 

E (<^m 2 L9»>o + \(g l m \Ag k m )v k + l(g l m \^ 2 \g k jv k 

k=l ^ 

N / 1 1 

E (9lW 2 \g k Jv k + ^J^AgtK + -M m \Ag k m )v k 



k=\ 



N 



= Y,(9 l m\V{»,v)\9 k J > 

k=l 
N 

= ^(g l m \v 2 v(^v)\g m ) , 

k=l 

N 

= Y,(9 l m\v 2 V{^v)\g m ) , 



(20) 



k=l 



where the index I — 1, . . . , N runs over all basis states 
and the notation g m = g m {y k ) is used. The potential 
V(/i, v) for the diamagnetic hydrogen atom is given in 
Eq. ©. All integrals in Eq. P0)) can be obtained an- 
alytically, and are presented in Appendix rAJ The set 
of equations (|20[) is a 3A-dimensional Hermitian positive 
semidefinite linear system for the coefficients v k , V k , V k , 



k = l,...,N, and must be solved, e.g., using a Cholesky 
decomposition [32| of the left-hand side matrix, at every 
time step when numerically integrating the equations of 
motion (|19[) . Technical remarks for the time propagation 
of coupled wave packets via the numerical integration of 
the Eqs. QU) will be given in Sec. InTCl 
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B. Hydrogen atom in crossed fields 



GWP (23) yields 



The rotational symmetry of the hydrogen atom in a 
magnetic field as discussed in Sec. IIII Al is broken when 
an additional electric field with a different orientation is 
applied. In crossed fields none of the three degrees of 
freedom can be separated. The paramagnetic term that 
contributed only a constant energy shift within the sub- 
space of constant m in the diamagnetic hydrogen atom 
must now be taken into account since l z is not con- 
served. The evolution of wave packets is determined 
by the time-dependent Schrodinger equation ([5]) with 
T given by minus one half times the Laplacian in the 
4D Kustaanheimo-Stiefel coordinates, and V defined via 
Eqs. (2) and (D as 



1 



V = au + -f3[(uip 2 - u 2 pi){u 3 + u 4 ) 
+ (u 3 P4 - u A p 3 )(u\ + u%)] 



(21) 



As trial functions for the time-dependent variational 
principle we use the superposition 



N 



^(z)=E<?(y & ) 



fc=i 



where 



g k ^g{y k ) = e* 



(22) 



(23) 



are the restricted Gaussian wave packets derived 
in the preceding paper (30l |. which depend on the 
5N time-dependent variational parameters y k = 
( r y k , a k ,a k ,a k ,a k ) (see Eq. (Jj), combined in the param- 
eter vector z = (y , . . . ,y ). The equations of motion 
for the variational parameters are obtained by evaluat- 
ing the TDVP in Eq. (T3) for the trial funct ion (22). 
The procedure is similar to that in Sec. IIII Al Letting 
the time derivative and the Laplacian act on a restricted 



+ itrA k )g k 



-uA k u~j k - 2u(A k ) 2 u 



2 uV 2" ) 9 K , 



(24) 



and defines a scalar Vq and a 4 x 4 matrix V k as the co- 
efficients of the polynomial in u for each GWP with k = 

k\2 



1. 



,N, i.e. 



itrA k -j k and V 2 k /2 = -A k -2(A k ) 



Since the special structure of the matrices A k in Eq. (7) 
is maintained in the squared matrices (A k ) 2 , that struc- 
ture carries over to the 4x4 complex symmetric matrices 
V 2 due to their definition in Eq. (24). Therefore, they 
have only four independent coefficients V k , V k , V k , and 
V k in the notation of Eq. (7). The equations of motion 
for the variational parameters y k = (7 , a k , a k 



, N can be written as 



A k 



-2{A 



k\2 



7 fe = i tr A k - v k 



1 
2 
■ 



V 2 k , 



k a k ) 

(25a) 
(25b) 



where the time-dependent parameters 

( v o > Kt! i Vxj Vy) are obtained at every time step by 
solving a linear set of equations. Using the derivatives 
of the restricted GWPs with respect to the variational 
parameters, 



19 



dg k 
da k 



i{u 2 + u 2 )g k , 



dg k ■/ 2 2^ k 9g k 
— = i{u 3 +u 4 )g , — j : = 2i{u 1 u 3 -u 2 UA)g 

OCL,, CJ CL „ 



dcf 

= 2i(uiu 4 + u 2 u 3 )g k , 

oa y 



(26) 



the required linear set of equations is derived from (15) 



J 



N / 1 1 

E\ jlk k 1 rife T/fe 1 jlk T/fe 1 jlkj/k i jlkirl 
I -'ll^O + 1 127^ V fi + J 137J V v + 1 14 V x + 1 15 V y 



k=l 
N 



k=l 
N 



k=l 



E ( J » v * + + ll ™l V » + ll & V * + ll £ V v 



E ( J i3«o + I l 2 k 3 \v k + I lk J-V k + I l 3 k V k + I^V X 

k=l ^ 

N ( 1 1 

E Il M + 0^ + T lh v » + i l M + i%y x 

k=l ^ 

N ( 1 1 

E + i&rf + i lk -v k + I lk V k + i£v t 



N 



= 

fe=l 

N 

— V T lk 

fe=l 
N 

— V T lk 

fc=l 
N 

- V T lk 
fe=l 

N 

- V T lk 



(27) 



fe=i 
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with I = 1, . . . , N. All integrals / in Eq. (gTJ) are denned 
and listed in Appendix [B] The potential (|22| for the 
hydrogen atom in crossed electric and magnetic fields, 
including the paramagnetic contribution, enters the inte- 
grals on the right-hand side of Eq. (|2"T|) . The linear set 
of equations ([27]) is Hermitian positive semidcfinite [see 
Eq. (|20p for the diamagnetic hydrogen atom] and can be 
solved using a Cholesky decomposition of the left-hand 
side matrix. 



C. Numerical time propagation of coupled wave 
packets 

An initial wave packet given as the superposition of 
basis states in Eq. (fl6|) or ([22)) can be easily propagated 
for the field-free hydrogen atom because the basis states 
remain uncoupled and the time-dependence of the ba- 
sis states is known analytically [30J. The external fields 
lead to couplings between the basis states, and the time- 
dependence of the variational parameters must be deter- 
mined numerically. The setup of the equations of motion 
has been discussed in Secs. lIII Al and lHTBl The numerical 
integration, however, of Eqs. (fT9f and (|25| is nontrivial 
and further remarks are necessary. 

1. Time propagation of the width matrices 

For better numerical performance it is advantageous 
pnl [33j to introduce, for each width matrix A, two aux- 
iliary time-dependent 4x4 matrices B and C in such a 
way that 

A = ^BC- 1 . (28) 

The equations of motion (|25ap and similarly Eqs. Q19ap 
and (|19bj) are then replaced with the equivalent differen- 
tial equations 

B k = -V k C k , 

C k = B k , (29) 

with the initial values B(0) = 2A(0) and (7(0) = 1. 
In the case of the diamagnetic hydrogen atom the ma- 
trices A and V2 are diagonal with diagonal elements 
{a p , a M , (V, a„} and {V^, V^, V u , V u }, respectively. The 
matrices B and C have the same structure, and thus 
the total number of parameters per basis state that 
must be integrated (including the scalar 7) increases 
from three parameters (7, a^,a„) to five parameters 
(j,b^,b„,c^,c v ). 

For crossed fields the increase of the number of param- 
eters is even more rapid. In that case the matrices B 
and C are no more complex symmetric. Without taking 
care of the special structure ([7]) of the matrix A the in- 
troduction of the B and C matrices would require the in- 
tegration of 32 complex parameters per GWP in the two 



matrices B and C instead of four complex parameters in 
the width matrix A. However, the special structure of 
the matrix A can be exploited to halve the number of 
independent parameters from 32 to 16 in the matrices B 
and C. Details are given in Appendix [Cj 

When integrating the equations of motion most of the 
computational effort is invested in solving the set of 3N 
linear equations (|2U)l or the 5N linear equations (|2T|) at 
each time step. The dimension of those equations is not 
affected by the introduction of the auxiliary matrices B 
and C, and thus the increase of the number of param- 
eters in the differential equations (|2"9"|) does not imply 
a significant increase of the total computing effort. In 
fact, due to the better numerical behavior of Eq. (|29p as 
compared to Eq. (|25ap and Eqs. (|19ap . (|19bp . larger step 
sizes of the numerical integration are possible and the 
total computing time is decreased. 



2. TDVP with constraints 

The equations of motion resulting from the TDVP es- 
pecially for a large number of coupled GWPs become 
badly behaved from time to time during the integration. 
In the general formulation of the TDVP at each time step 
the linear set of equations (fT3"|) must be solved for the 
equations of motion of the variational parameters, i.e., 
the time derivatives z. In the course of integration, de- 
pending on the number of coupled GWPs, it will happen 
sooner or later that the matrix K in Eq. (|13p associated 
with the set of linear equations becomes ill-conditioned, 
or even numerically singular. As a result the time step of 
the integration routine becomes extremely small, render- 
ing the method of GWP propagation impracticably slow. 
In the worst case the wave packet propagation can stick 
completely. 

Matrix singularity problems arise from overcrowding 
the basis set, i.e., from situations where fewer GWPs 
would be sufficient to represent the wave function. On 
the other hand for an accurate approximation of the wave 
function it is desirable to have a large number of ad- 
justable parameters. However, there is a discrepancy be- 
tween the number of GWPs necessary to yield accurate 
results and the maximum number of GWPs that can be 
propagated using the TDVP without numerical difficul- 
ties [3J| . There exist different proposals to overcome this 
numerical problem [HI [H, H, EH, |37|, H IM S3 • Here we 
adopt the constrained time-dependent variational princi- 
ple [401 ] , where inequality constraints of the form 

/fc(z,Z*) EE f k (z r ,Zi) EE /fe(z) > fk,mm , 

f k e R, k = 1,2,3,... (30) 

are taken into account in the variational process, and 
complex quantities are split into their real and imaginary 
parts, denoted by the subscripts r and i, respectively, and 
thus z ee (z r ,Zj). The functions fk must be chosen in 
such a way to prevent the matrix K from becoming sin- 
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gular. As long as /fc(z r , > fk,min f° r au au param- 
eters evolve according to Eq. (fT3"|) without being affected 
by the constraints. However, when /fc(z r ,Zi) = /fe )m i n 
and /fe(z r , Zj) < for, say, fc = 1, . . . , j we introduce La- 
grangian multipliers and obtain an extended set of linear 
equations 



K 


M T 


M 






""■*-(£; t)' fi -u)- (31) 

where the matrix if and the vector h are the complex 
quantities of Eq. (| 13(1 . The Lagrangian multipliers are 
A € MP and M = || with f = (/i, . . . , fj) is a real valued 
j x 2n p matrix. Details of the implementation of the 
constrained TDVP are given in 40]. If no constraint is 
active, i.e., j = 0, then Eq. (f3Tj) obviously reduces to the 
real formulation of Eq. fT3|) . 

IV. RESULTS AND DISCUSSION 

In this section we present examples for the fictitious 
time wave packet propagation of the hydrogen atom in 
external fields. Autocorrelation functions between the 
initial and time propagated wave packets are computed. 
Quantum spectra are obtained by the frequency analy- 
sis of the autocorrelation function and compared with 
numerically exact diagonalizations of the Hamiltonian. 

It turns out that a sensible choice of an appropriate 
initial state if>(Q) is crucial for the successful application 
of the TDVP. For an unreasonable choice the numeri- 
cal problems discussed in Sec. MI C 21 occur for few ba- 
sis states already, and bad, unconverged results are ob- 
tained. The conventional way to construct an initial wave 
packet by placing a certain number of unrestricted GWPs 
at various positions in coordinate and momentum space 
is not possible for the restricted GWPs in the KS coor- 
dinates. In the calculations of the diamagnetic and the 
crossed fields hydrogen atom we achieved optimal results 
by first choosing only one 2D or 3D Gaussian wave packet 
in the physical coordinates, which was then expanded in 
a set of N restricted GWPs as explained in Ref. [3(| for 
the field-free hydrogen atom. The external fields lead to 
couplings between the basis states and imply a compli- 
cated time development of the initial state as compared 
to the field-free hydrogen atonijWhere the wave packet 
propagation is periodic in time |30l | . 

A. Diamagnetic hydrogen atom 

The initial wave function is most conveniently chosen 
to be a GWP in parabolic coordinates 

Tj) = ^ e -(«-«o) 2 /(4<r 2 )-(»?-»?o) 2 /(4^ 2 )+ip ! : (C-Co)+ip™(';- 

(32) 




FIG. 2: (Color online) Fictitious time evolution of the state 
(]32|) with po — 6.0, zq = and a nonzero initial mean mo- 
mentum. The wave function is plotted for different values of 
the dimensionless fictitious time t. The initial wave packet 
gradually becomes delocalized. Lengths are given in scaled 
atomic units n c sao with ao the Bohr radius (see Eq. @). 



with center (£o>??o); width a and mean momentum 
(P(oiPvo)- The GWP is expanded in terms of the ba- 
sis states (fTT)) according to the procedure described in 
detail in Ref. [30J, including the Monte Carlo technique 
with importance sampling. The procedure yields the ini- 
tial values of the variational parameters 7 fe ,ajj,a^, k = 
1,...,N. 

However, it is not realistic to propagate several thou- 
sands of basis states numerically with the full coupling. 
Reliable results are obtained by far fewer basis states 
than used in the expansion and propagation of the GWP 
(f32|) in the field- free hydrogen atom [jjjj. Reasonable 
numbers of basis states are in the range of N=10-100. A 
,numerical example is presented for the magnetic quan- 
tum number m = 0, where N = 70 basis states are used 
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20 40 60 80 100 



x 

FIG. 3: (Color online) Real part of the autocorrelation func- 
tion (^(t) = <V o t (0)|V'6 t (' r )> for the GWP <E2} with the cen- 
ter po = 6.0, zo = 0. (a) Signal of the projected state with 
even parity, and (b) odd parity. The fictitious time r and the 
signal C(t) are in dimensionless units. 



for the expansion and propagation. The damping factor 
e is set to e = 0.1 

Each basis state has three variational parameters 
7 fc ,a^,aj5, and therefore N basis states require the so- 
lution of a 3iV x 3iV matrix equation after every integra- 
tion step, and the usual numerical problems mentioned in 
Sec. IIIII occur with increasing number of basis states. It 
turns out that constraints on the imaginary parts of the 
phase parameters of the form Im.7 > j m i n — —4.5; k 
1 Y are suitable to regularize the equations of mo- 
tion with regard to a fast integration. These constraints 
present simple lower bounds on the amplitudes of the 
wave packets and avoid matrix singularities caused by 
extremely large overlapping wave packets. 

The accuracy of the expansion ([3"2"| of the GWP with 
only N = 70 basis states is very good. The time evo- 
lution of the wave function is shown in Fig. [21 The 
probability density p\ip(p, z)\ 2 for six different times r = 
0.4, 0.8, 1.2, 3.0, 5.0, 7.0 is shown. The parameters of the 
potential in the Hamiltonian ([9]) are set to a = 0.5 and 
f3 = 0.2. The tt periodicity of the evolution of the wave 
function that is present in the field-free hydrogen atom, 
is destroyed now. 

The autocorrelation function of the propagation can 
be used to extract spectral information by Fourier trans- 
formation or harmonic inversion jJJ, H2, IH, 0, [4|| of 



(a) 0.3 r 
-a 0.2 - 



T 0.2 - 

I " - 

* - 



2 4 6 8 10 12 14 
n eff 

FIG. 4: (Color online) Spectra with (a) even and (b) odd z- 
parity extracted from the autocorrelation function C ± (r) = 
(ipo(r = 0)\ipQ (r)) computed from the evolution of the wave 
function (|32|) plotted in Fig. [2] The amplitudes are given by 
the magnitude of overlap between the initial wave function 
and the respective eigenstates. For comparison the positions 
of the numerically exact eigenvalues obtained from a diago- 
nalization are plotted with blue lines in the lower panels of 
the figures. The related eigenenergies and the magnetic field 
strength follow simply from Eq. Q. The effective quantum 
number n e ff and the overlap matrix elements are in dimen- 
sionless units. 



the time signal. The center of the Gaussian (|3"2")l is 
po = 6, zo = and the initial mean momentum is chosen 
in such a way that states around an effective quantum 
number of n e g « 6 are excited. 

To reduce the density of states the autocorrelation 
function is separately computed for the subspaces of even 
and odd parity by taking the symmetrized and antisym- 
metrized states ip^{p, z) = ipo(p, z) ± ipo{p- — z). The au- 
tocorrelation function C ± (r) = ('0o = (O)|'0o = ( T ))7 is shown 
in Fig. [3ja) for symmetrized states and in Fig. [3^b) for 
the antisymmetric states. The spectral results for the 
diamagnetic hydrogen atom, obtained from the time sig- 
nals are plotted in Fig. 2] A harmonic inversion has been 
employed. The amplitudes of the peaks are determined 
by the magnitude of the overlap between the eigenstates, 
denoted by |n e jf), and the initial states ^(O) in Fig.0Ja) 
and Fig. [31(b) , respectively. The amplitudes are plotted 
with red lines. The numerically exact eigenvalues of the 
diamagnetic hydrogen atom are plotted with blue lines 
for comparison. The agreement of the positions is excel- 
lent. The highest amplitudes are located in the region 
n e g ~ 6 according to the choice of the input parame- 
ters of the initial GWP in Eq. . The multiplicity of 
the states with even or odd z-parity resulting from the 
same principle quantum number n is determined by the 
number of positive and negative eigenvalues (— \) l+m of 
the z-parity operator acting on the spherical harmonics 
Yi m (6, ip) with I < n. 

The values of the parameters a = 0.5 and j3 = 0.2 used 
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FIG. 5: (Color online) Effective quantum numbers n e g at the 
field-free ionization threshold E = a = for states of (a) even 
and (b) odd parity. The propagation involves N = 90 basis 
states. The variational results (red lines in the upper panels) 
are in excellent agreement with the exact time-independent 
results (blue lines in the lower panels). The effective quan- 
tum number n e ff and the overlap matrix elements are in di- 
mensionless units. 



for this computation still present a mainly harmonic sys- 
tem with a perturbation for low energies. As mentioned 
above the dynamics of wave packets is exact to all or- 
ders in the field strengths within the allowed set of trial 
wave functions, i.e., the variational approximation only 
concerns the restriction of the Hilbert space. Therefore, 
the method is not restricted to the pertubative regime 
but even allows for the computation of eigenvalues in the 
strong anharmonic regime. Fig. [5] presents results at the 
field-free ionization energy E = a = and j3 — 0.5 for (a) 
even parity and (b) odd parity. A number of N = 90 ba- 
sis states was used in the computation. In the presence of 
the magnetic field these states at the field-free ionization 
energy E = are still bound. The agreement between 
the eigenvalues computed variationally (red lines) and 
the numerically exact results (blue lines) is very good. 
The related field strengths are easily obtained from Eq. 
(U]) by B — p/rv^Q. The underlying initial wave packet 
(f3"2| is initially centered at po = 4.39, Zq = 1 and has 
zero mean momentum. As mentioned above the posi- 
tion, momentum, and width of the initial GWP deter- 
mine the spectral region for n e g where strong peaks are 
expected. However, within that region some eigenstates 
\n e g) can be near orthogonal to the initial GWP and thus 
have nearly zero amplitude. Indeed, some lines are lack- 
ing in the variational computation. The missing states 
can be revealed by choosing several initial GWPs, which 
have larger overlap with those states. An example of the 
influence of the chosen initial GWP on the peak ampli- 
tudes will be given in Scc. lIVBl for the hydrogen atom in 
crossed electric and magnetic fields. 
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FIG. 6: (Color online) Spectra with (a) even and (b) odd 
z parity of the Hamiltonian ([2| with a = 0.5, /3 = 0.05, 
£ = 0.01 obtained from the propagation of two different 3D 
GWPs. Green and red line (upper panels in the figures): 
x = (6,0,0), po = (0,±1/V5,1/V2), respectively. The 
eigenvalues are extracted from the autocorrelation function 
by Fourier transformation. The peak positions agree very 
well with the numerically exact eigenvalues of the effective 
quantum number marked by blue lines in the lower panels of 
the figures. The related eigenenergies and the field strengths 
follow from Eq. Q. The effective quantum number n e s and 
amplitudes are in dimensionless units. 



B. Hydrogen atom in crossed fields 

For the hydrogen atom in crossed electric and magnetic 
fields the propagation of 3D GWPs is computed starting 
from the time-dependent Schrodinger equation ([2]) with 
parameters a — 0.5, — 0.05, and C = 0.01 in Eq. (g]). 
The choice of an appropriate initial state "0(0) is very im- 
portant for the successful application of the TDVP. We 
achieved optimal results by first choosing one 3D Gaus- 
sian wave packet in physical Cartesian coordinates with 
the center Xo and width a in position space and center 
Po in momentum space 



^(x) = (27TCT 2 ) 



-3/4 



exp 



x ) 2 



4cr 2 



+ ipo • (x - x ; 



(33) 

which is then expanded in a set of N restricted GWPs. 
The external fields lead to couplings between the basis 
states, and the time-dependence of the variational pa- 
rameters must be determined by the numerical integra- 
tion of Eq. (|25p . For better numerical performance we 
resort to the TDVP with constraints [401 ] mentioned in 
Sec. IIII C 21 As for the diamagnetic hydrogen atom con- 
straints of the form Ini7 fc > 7 m ; n = —4.0, k = 1, . . . , N 
are imposed on the imaginary parts of the phase parame- 
ters j k . Once a time-dependent wave packet ([6]) is deter- 
mined the eigenvalues n e ff of the stationary Schrodinger 
equation @ are obtained by the frequency analysis of 
the time signal (|10p with the amplitudes Cj depending 
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on the choice of the initial wave packet. In perpendicular 
crossed fields the z parity is conserved. Spectra with even 
and odd z parity obtained from the Fourier transforms 
of the autocorrelation functions C ± (t) = ( , ± (O)|?/' ± (t)) 
of the parity projected wave packets are shown in Fig. 
[6] In Fig. [6]Ja) the eigenvalues with even parity and in 
Fig. Hlb) the eigenvalues with odd parity are plotted. 
The green and red lines result from the propagation of 
two different 3D GWPs with a = 3.5, e = 0.15 and the 
same initial position x = (6,0,0) but different initial 
mean momenta po = (0, ±l/\/2, l/\/2), respectively. A 
number of N = 41 and N — 31 basis states were coupled 
in the calculations. The line widths, i.e., the resolution 
of the spectra, is determined by the length of the time 
signal r max . The eigenvalues obtained by numerically 
exact diagonalizations of the stationary Hamiltonian (|2|) 
are shown by the blue lines. The line- by- line compari- 
son shows good agreement between the exact spectrum 
and the results obtained by the wave packet propagation. 
The amplitudes of levels indicate the excitation strengths 
of states with higher or lower angular momentum l z by 
the two initial wave packets rotating clockwise or anti- 
clockwise around the z axis. 



V. CONCLUSION 

The Gaussian wave packet method is known to be well 
suited for systems with nonsingular smooth potentials 
but not so for systems with singular potentials such as the 
Coulomb potential. Therefore so far it failed when ap- 
plied to atomic systems. Using the Kustaanheimo-Stiefel 
regularization of the Coulomb potential and introducing 
a fictitious time variable we have now made it applica- 
ble to atomic systems by using restricted GWPS and the 
time-dependent variational principle. The special appeal 
of the GWP method lies in the fact that relatively low 
numbers of time-dependent basis states are sufficient to 
derive the spectrum as compared to time-independent 
matrix diagonalizations. The advantage of using the fic- 
titious time is that the computations are exact and ana- 
lytical for the field-free hydrogen atom, which means that 
for perturbed atomic systems only the deviation of the 
potential from the Coulomb part must be taken into ac- 
count in the variational approximation. We have shown 
that the method can be especially adapted for systems 
with, e.g., cylindrical or spherical symmetries. 

Quantum spectra of the hydrogen atom in static ex- 
ternal fields can nowadays be computed quite efficiently 
by matrix diagonalization of the Hamiltonian in a suffi- 
ciently large basis set, and thus the method proposed in 
this paper might appear to be rather specific as an al- 
ternative tool for studying this system with complex dy- 
namics. However, quantum computations for many-body 
Coulomb systems are certainly a nontrivial task. The 
topic of wave packet dynamics in systems with Coulomb 



interactions covers a large body of problems ranging from 
atomic physics to physics of solid state, where Coulomb 
interaction plays an important, often crucial, role. In 
many-body physics, in particular, in the physics of solid 
state, theoretical methods well suited for studying the ef- 
fects stemming from Coulomb interactions are still lack- 
ing. The majority of the available methods, e.g., the 
method of pseudopotentials in atomic physics and the 
Fermi and the Luttinger liquid theories for solid conduc- 
tors, are basically indirect and substantiated neither from 
the theoretical nor from the experimental side. For this 
reason they still remain, to a certain extent, disputable. 
In the present paper we have successfully applied the 
Gaussian wave packet method to Coulomb systems with 
two and three nonseparable degrees of freedom. If the 
method can be further extended to larger systems with 
more degrees of freedom it will allow for a wide range of 
future applications in different branches of physics. 

APPENDIX A: INTEGRALS FOR THE 
DIAMAGNETIC HYDROGEN ATOM 

With the basis functions g m defined in Eq. (flT]) , and 
using the notation = a* — (a 1 )*, a v = a£ — (a l v )*, 
7 = 7 fc — (7')*, and m > 0, the absolute value of the 
magnetic quantum number the integrals in Eq. (|20[) take 
the form 

(sll/MISm) 

POO poo 

= Att 2 dfi d i /(H 2m+1 /(M 2 ^ 2 )e''^ 2+0 " I/2+7 ) , 
Jo Jo 

(Al) 

where f(p 2 ,v 2 ) is a poynomial in a 2 and v 2 . The in- 
tegrals can be factorized, and with x — u 2 or x — 
v 2 the products basically take the elementary form 
/ °° x n e~ ax dx — -£tt for integers n > and Re a > 0. 
The integrals on the left-hand side of Eq. (|20|) read 



(din 




n 2 (m\) 2 ^ 


— c 


9rn) = 


(-a f ,a u ) m + 1 




9rn) = 


-4— (m + 1) , 




{glW 2 


<?m) = 


— — (m + 1) , 

— %a v 




(sLIm 4 


9 m ) = 


c 

-(m + l)(m 

-< 


+■ 2) 




9 m ) = 


C (m+1) 2 




<5> 4 


k \ 
9,n) = 


Q 

-(in + l)(m - 

-at 


+-2) 



With the potential V(/i, v) given in Eq. the integrals 
on the right-hand side of Eq. PD|) are obtained as 
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{9 l m \V(fx,v)\gi) = -^( atl + a,)(l + m)[(2 + '3m + m 2 )(3 2 + 8a^a} , 
(g l m \n 2 V(fi, ) = — J-j (m + 1){(1 + m)(2 + m)[o p (2 + m) + a„(3 + m)]/3 2 + 8o M o v [o^ + 2a,, + (a M + a u )m]a} , 
(g l m \v 2 V([i, ) = -^-3- (m + 1){(1 + m)(2 + m)[a p (3 + m) + a„(2 + m)}/3 2 + 8o p o„[2a M + + (a M + a„)m]a} . 

APPENDIX B: INTEGRALS FOR THE HYDROGEN ATOM IN CROSSED FIELDS 

The integrals in the linear set of equations (|27[) take the form 

<fl'|/(u,V u )|fl*> 

= J d 4 ue- l i u{Alru+ ^f(u,V u )e^ uAku+ ^ . (Bl) 
With the notation A = A k — (A 1 )*, 7 = j k — (7')* the integrals on the left-hand side of Eq. (|27| simplify to 

4" = (ff'l/i/ils*) = / d^ufjje^^ , (B2) 



with /1 = 1, /2 = u\ + 1*2, /3 = u§ + u|, fi = U1U3 — U2U4, and — uiu^ + 1*21*3. The integrals have the properties 
J|* = and 7^ fe = Using c = 7r 2 e^ and /i = 1/V- detA = l/(a 2 + a 2 - a^a,,) we obtain 



7 X1 = ft,c , I l 12 = —ia v h 2 c , 7 13 = —ia^h^c , 
/{I = 2i fl2; /i 2 c , 7(| = 2ia 9 /i 2 c , J^j = -2a 2 u h 3 c , 
4s = ~(a^ + a * + Oy)^ 3 c , 7m = 4a„a x h 3 c , 
7^ = ±a v a y h 3 c , = -2a 2 /i 3 c , = i afl a x h 3 c , 

r' fc _ /f„ „ f,3„ jlk _ n/2 q„2 „ „ u3. 



7^ = -8a x a y h 3 c , I l 5 k 5 = 2(a 2 - 3a 2 - a,,a u )h 3 c . (B3) 



7^ = 4a M a y /i 3 c , 7^ = 2(a 2 - 3< - a^aJjhTc , 

The integrals on the right-hand side of Eq. (|27|) are defined as 

45 = (ffU^ls*} ■ (B4) 

The potential V defined via Eq. ^ can be split into its harmonic and diamagnetic part, 

V a = au 2 + \p 2 {u\ + u 2 2 ){ul + u 2 )u 2 , (B5) 

o 

and the terms of the paramagnetic and electric field contributions, 

V b = ^(3[{u lP2 -u 2Pl ){u 2 + ul) (B6) 

+ ("3P4 - U 4 p 3 )(ul + ulj] + C(UXU 3 - U 2 Ui)u 2 . 

Note that the paramagnetic term in Eq. (|B7[) contains derivatives with respect to the KS coordinates and thus the 
integrals must be solved by application of Eq. (|B1|) ■ We obtain 



4la = i(a^ + a v )[(a f ,a„ + 2(al+a 2 y ))/3 2 h 2 /A + a}h 2 c 1 

4aa = i 2a l a l + ( a l + a l)i Qa t + 2 ( a l + a l)) + a»av{3al + 8(a 2 + a 2 y ))]p 2 h 5 c/4 + [a^a^ + 2a v ) + a 2 + a 2 y ]ah 3 c , 

4,3a = [3a^a„ + 8a^a u (a 2 x + a 2 ) + 2(a 2 + a 2 ) 2 + a 2 (2a 2 + 9(a 2 + a 2 ))]/3 2 /i 5 c/4 + [a^(2a M + a v ) + a 2 + a 2 y ]ah 3 c , 

4*0 = ~( a n + a y )a :c [3(a A1 a^ + a 2 +a 2 )/3 2 /i 2 +4a]/i 3 c , 

4L = -{a^ + a„)a y [3(a^ + a 2 x +a 2 y )p 2 h 2 + 4a]/i 3 c, (B7) 
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1% = 2(a M + a„)(a k y a x (3 - a k a y f3 + a x ()h 3 c , 
I-v2b = — 2i(2a M a,/ + 3a^ + + a 2 )(a k a x (3 - a k a y (3 + a x Qh A c , 
1% = ~2i{3al + 2a^a„ + al + al)(a k y a x p-a k x a y f3 + a x ()h i c, 
I lk Ab = 2i(a^ + a v ) [-6a k a x a y /3 + a^a^ + ba x - a 2 )(i + (a^ + 5a 2 x - a 2 )(]h 4 c , 
I l v % b = -2i(a /t + a w )[a*(a |i a I ,- al + 5a 2 y )f3 - 6a x a y (a k f3 + (^c , (B8) 

I 



The right-hand side vector in Eq. 
corresponding terms in Eqs. jB7| 

jlk 

vja 



27)) is the sum of two 
and (JBSj), i.e., I lk = 



APPENDIX C: STRUCTURE OF THE 
MATRICES B AND C 

A structure of the matrices B k and C k is searched 
which is preserved in the matrix product V k C k in Eq. 
(|2"5)) . where V k has the same structure as A in Eq. ([7]). 
This is provided by the form 



B = 



11 


&12 


&13 


&14 \ 


12 


&11 


&14 


-&13 


31 


&32 


&33 


&34 


32 


-&31 


-&34 


^33 / 



\ b 3 2 - 



C 
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Cll 


C12 


Cl3 


Cl4 ^ 




-C12 


Cll 


Cl4 


-C13 




C31 


C32 


C33 


C34 


\ 


C32 


-C31 


-C34 


C33 / 



(CI) 



as can easily be shown by explicit multiplication. The 
superscript k running over all GWPs has been omitted 
here. Compared to Eq. (JT]) the number of independent 
parameters per matrix increases from 4 to 8, however, 
this is still less than 16 parameters for a general 4x4 
matrix without any special structure. 

The matrix A = \BC~ X can be calculated analytically. 
To this end we introduce the auxiliary matrix 



D 



V 



C33 
C34 
-C31 

-C32 



-C34 
C33 

-C32 
C31 



-C13 -C14 \ 
-C14 C13 

Cll -C12 
Cl2 Cn / 



(C2) 



The product C x = CD yields 



C x = 



( h k 

-k h 

h 

\ k 




-k 
hj 



(C3) 



with 



-C13C31 - C14C32 
C14C31 - C13C32 



C11C33 + C12C34 , 
C12C33 - C11C34 . 



The matrix C\ can be easily inverted, and thus allows for 
the calculation of A = \BC~ X = \BDC^ X . With h' = 
h/[2(h 2 + k 2 )}, k' = k/[2{h 2 + k 2 )] the four independent 
parameters in Eq. ([7]) read 

<V = (&IIC33 - &14C32 + ^12C 3 4 - &13C3l)^' 

+ (6i4C 3 i + 612C33 - 611C34 - 613032)^' , 

a v = (>33Cll + ^34Cl2 - &31C13 - b 32 C U )ti 

+ (&33Cl2 - &34C11 - 632C13 + 63lCl4)fe' , 

a x = (&13C11 + 614C12 - &11C13 - h2Cu)h' 

+ (&13Cl2 - &14C11 - &12C13 + buCi4)k' , 
a y = (614C11 - 613C12 + &12C13 - buciijh' 

+(&i3Cn + 614C12 - 611C13 - b 12 ci 4 )k' . (C4) 
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